The decay and collisions of dark solitons in superfluid Fermi gases 
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We study soliton collisions and the decay of solitons into sound in superfluid Fermi gases across 
the Bose-Einstein condensate to Bardeen-Cooper-Schrieffer (BEC-BCS) crossover by performing 
numerical simulations of the time-dependent Bogoliubov-de Gennes equations. This decay process 
occurs when the solitons are accelerated to the bulk pair-breaking speed by an external potential. A 
similar decay process may occur when solitons are accelerated by an inelastic collision with another 
soliton. We find that soliton collisions become increasingly inelastic as we move from the BEC to 
BCS regimes, and the excess energy is converted into sound. We interpret this effect as being due 
to evolution of Andreev bound states localized within the soliton. 



I. INTRODUCTION 



Since the realization of Bose-Einstein condensation in dilute gases, experimental and theoretical research has shown 
that solitons play a key role in their dynamics (see Ref. [i| for recent reviews). Solitons have been generated in a 
wide variety of contexts such as, for example, phase-imprinting in elongated Bose-Einstein condensates (BECs) 0, 
0|, transport in optical lattices 0, transport past single defects or through disordered potentials BEC 
collisions [9l-[Tl| and interferometry [l^ . Consequently, solitons might be expected to play an equally important 
role in the dynamics of degenerate Fermi gases. Recent work has predicted the existence of black 2M ^^'^ S^^Y 
solitons jl5l - [l7| across the Bose-Einstein condensate to Bardeen-Cooper-SchriefFer (BEC-BCS) crossover, and that 
they may perform stable oscillations in a trap [l^. However, if solitons in Fermi gases are to be as ubiquitous as 
solitons in Bose gases, they must be robust in the presence of defects, rapidly-varying potentials, superfluid flow, 
sound and other solitons. Can they be easily produced and are they long-lived? If not, what excitations will be 
observed? 

In this paper, we address some of these questions by solving the time-dependent Bogoliubov-de Gennes (TDBdG) 
equations across the BEC-BCS crossover. We confirm a previous hypothesis that soliton solutions do not exist with 
a speed v above the bulk pair-breaking speed Vpb [l6| . This means that on the BCS side of unitarity, where Vpt, is 
smaller than the sound speed c, the maximum v is Vpb, in accordance with the Landau criterion. In contrast, in the 
BEC regime the maximum v is c, as predicted by the Gross-Pitaevskii equation If a soliton in the BCS regime 

is accelerated to Vpb by an external potential, it abruptly disappears and its energy is converted into sound. We stress 
that this decay process is distinct from the snake instability [20-2^, by which a soliton can decay into quantized 
vortices in a three-dimensional superfluid. In our calculations, in fact, we assume that the dynamics are restricted 
to longitudinal motion only, the gas remaining uniform in the transverse directions. This is equivalent to assuming 
that, in a real experiment, the superfluid is sufficiently tightly-confined in the transverse directions to suppress the 
snake instability. Moreover, the predicted decay process occurs at zero temperature and is distinct from other decay 
mechanisms considered previously (24i-i26i] . 

We analyse this decay process by looking at the low-lying eigenstates of the Bogoliubov quasi-particle spectrum. 
The lowest state is found to be localized in the vicinity of the soliton, with an energy below the gap for the continuous 
spectrum of the extended Bogoliubov states. The origin of this localized state (Andreev state ,27n29i] ') is the fact that 
the energy cost for creating a fermionic excitation near a minimum of the order parameter is reduced with res pec t 
to the bulk value. The Andreev states associated with dark solitons at rest have already been discussed in Ref. jl4|. 
Here we find that Andreev states also exist for solitons with finite v. However, as v approaches Vpb, the excitation gap 
in the continuous spectrum approaches zero and the Andreev bound state is lost. We also study the soliton energy 
Es as a function of v (the energy dispersion), and find that the soliton still has a finite Eg and phase jump across 
it at the maximum v. The energy dispersion shows interesting structure as v approaches Vpb, in particular a local 
minimum. We explain that, in a real experiment, this local minimum would cause the soliton to decay at the value 
of V in the minimum rather than at the maximum v where the energy dispersion truncates. Moreover, the presence 
of a harmonic trap causes a further reduction in the value of v at which the soliton decays. This finite-size effect is 
reduced as we weaken the trap. 

The existence of a maximum soliton speed, or minimum soliton energy, has important consequences for the physics 
of soliton collisions across the crossover. In the BEC limit, due to the integrability of the one-dimensional Gross- 
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Pitaevskii equation for a homogeneous gas, sohton colhsions are elastic and hence solitons always emerge from a 
collision with the same velocity and energy at which they collided [s^, HH- However, we find that soliton collisions 
become increasingly inelastic as we reduce l/kja. As for the decay of a single soliton, fermionic quasi- particles are 
found to play a key role also in these inelastic collisions. This is consistent with the fact that the same collisions between 
solitons are found to be almost completely elastic |32j when described by means of a nonlinear Schrodinger equation, 
which is often used to investigate the dynamics of superfluid fermions in the BEC-BCS crossover, but includ ing only 
bosonic degrees of freedom. Since, counter-intuitively, fast solitons have less energy than slow solitons [l^, [3 Il9j |. 
the inelastic collisions cause the solitons to accelerate. The energy lost by the solitons is dissipated as sound. We 
interpret this effect as being due to the evolution of the Andreev bound states as the soliton changes shape during the 
collision; this evolution eventually causes an energy transfer from the soliton to the continuum of bulk excitations. If 
the solitons lose so much energy that their energy after the collision would be less than the minimum soliton energy, 
they are destroyed by the collision. 



II. METHODOLOGY 



We consider a three-dimensional superfluid Fermi gas with equal populations of the two spin components. We 
model its dynamics across the BEC-BCS crossover by solving the TDBdG equations [l^, HI] 

'H A(r,t) 
_A*(r,t) -H 

where H = — /i^V^/2m + U — m(0), in which m is the atomic mass, U is the external potential and /i is the 
chemical potential. The order parameter is calculated as A(r, = ~5 X^jj "'^'j'^j?! i'^ which g is given by ^/kfa = 
SnEf /{gkj:) + y^iEc/ (n'^Ef) [3^]. Here a is the 3D s-wave scattering length characterizing the interaction between 
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atoms of different spins, while Ef = ?i^fcj/2m and kf = (Stt^tt.) are the Fermi energy and momentum of an ideal 
Fermi gas of density n, respectively. The cut-off energy Ec is introduced in order to remove the ultraviolet divergences 
in the BdG equations with contact potentials. The density of the gas is n{r, t) = 2 lwj)(r, t)\^ . We impose that the 
potential U has no y or z dependence, and hence we may write the functions u^(r, t) and w,,(r, t) as u^(a;, t)e*('^!'^+'^^^) 
and Vri{x,t)e'^'^^^y'^^'^^ respectively, in which ky and k^ are quantized according to ky = 2'Kay/L±^ and k^ = ^-KazjL^^ 
where Uy and are integers and i jis the width of the box in the y- and z-directions. As initial states at t = 0, we 
find stationary solutions of Eq. ([Ij [13, [IB] • 

We also search for solutions of Eq. ([T]) satisfying A (x, t) — A[x — vt), given that ?7 = 0, by solving the equation 

' H^ + thvf^ A(0 

where ^ = x — vt, = —h^ /2m [9^/9^^ ~ ~ ^z] ~ M and is the energy of level rj. We use a generahzed secant 
(Broyden's) method to find self-consistent solutions [la]- Hence we find travelling soliton solutions in a homogeneous 
gas which are stationary in the frame of the soliton, thus generalising to the BEC-BCS crossover the solutions for 
solitons propagating in a homogeneous Bose gas O, [l^, [l^- This technique enables us to determine the soliton 
properties more accurately in the homogeneous gas and identify effects due to the trapping potential by comparison 
with the time-dependent simulations. 
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III. SOLITON DECAY 

Some of us have shown in a previous publication (l6| that the energy dispersion of the soliton in the BCS regime 
truncates below the speed of sound c. This result suggested that, more generally, soliton solutions do not exist 
for speeds greater than Vpb or c, whichever is the smaller. To support this hypothesis, and to show how the effect 
would manifest itself in a real experiment, we now present time-dependent simulations of the Bogoliubov-de Gennes 
equations [l^ [s^ in which we accelerate the soliton above the critical velocity. 

Firstly, for comparison, in Fig. [T] we present a simulation of a stable soliton oscillation in a trap at 1/fc/a = —0.5. 
As in a previous publication [l5|, we consider a soliton oscillating in a '^'^K superfluid, contained in the harmonic 
trapping potential U{x) = with uJx = 2tt x 50 rad s~^, L± = 3.3 /zm and a peak density np = 1.8 x 10^^ 

m^-^. Figures [TJa) and (b) show the density profile and phase of the order parameter. Note that only the region of 
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FIG. 1: (a) Grey-scale plot of density profile n{x,t) (black=high) and (b) phase of the order parameter A{x,t) for a stable 
soliton oscillation at 1/fc/a = —0.5 in a '^'^K superfluid, with ujj, — 2tv x 50 rad s~^, L± — 3.3 fim and a peak density 
Up = 1.8 X 10^* m"''. The soliton begins at rest at a distance Xo = 3.3 /im from the trap centre, (c) & (d): Corresponding 
plots for a larger Xo of 5.4 fim producing soliton decay. 

cloud near the center of the trap is shown, the low density tails of the cloud are outside of the field-of-view. The 
soliton begins at rest at a distance Xq = 3.3 fim from the trap centre, and is consequently accelerated by the harmonic 
potential. As this happens, the density profile of the soliton becomes shallower and the phase jump across the soliton 
reduces from tt. However, as the soliton passes the center of the trap and reaches its maximum speed, it is still 
localized in the density profile and the phase jump is well-defined. The soliton then begins to decelerate as it climbs 
the trap potential, and the phase jump increases towards tt. If the simulation were allowed to continue the soliton 
would perform a complete oscillation in the harmonic potential. 

Figures [TJc) and (d) show the corresponding density profile and phase for the larger Xq of 5.4 fim. The soliton 
begins to accelerate, and initially the soliton remains localized in the density profile with a clear phase jump. However, 
when the soliton reaches a critical speed it rapidly spreads out in the density profile and the phase jump disappears. 
At the end of the simulation we see only low-amplitude modulations of the density profile. By measuring their speed 
we identify these modulations as sound. 

To gain a deeper physical understanding of this effect, we now study the soliton in the homogeneous gas using the 
BdG equations ([2]). The magnitude of the order parameter |A(a;)| is plotted in Fig. [2{a) for 1/kfa — —0.5 and v = 
(sohd blue curve), 0.2c (solid green curve with circles), 0.3c (dashed red curve) and 0.38c (dot-dashed black curve). 
As V approaches the pair-breaking velocity Vpb = 0.41c, the soliton becomes very shallow and broad. This is reflected 
in the profile of the lowest Andreev state, shown in Fig. [2fb). At v = 0, |vo(a;)|^ is tightly localised at the soliton 
position (blue solid curve). Also note that there are many nodes in the function. As v increases to 0.2c (green 
solid curve with circles) and 0.3c (red dashed curve), the function widens and the nodes becomes gentle oscillations. 
When V reaches 0.38c, the function is broad and smooth with almost no discernable oscillations (black dot-dashed 
curve) . 

As V approaches the pair-breaking velocity Vpb = 0.41c for l/kfU = —0.5, the excitation gap in the continuous 
spectrum of the extended Bogoliubov states, shown by the black solid line in the Fig. [2]inset, approaches zero. Below 
the continuous spectrum, there is the single Andreev bound state, which is plotted for various v in Fig. [5Jb). As the 
excitation gap closes, the lower boundary of the continuous spectrum approaches the energy of the Andreev bound 
state, shown by the dashed line in the Fig. [5] inset. Note that the energy of the Andreev bound state is almost 
independent of velocity. Eventually, the localized Andreev bound state is lost and becomes an extended state. 

The loss of the Andreev state can be understood by a simple analytic argument. We treat Eq. ([2]) semiclassically by 
considering the situation where the order parameter A varies slowly over scales of the order of kj^ p7| . This is true 
for the soliton near its maximum velocity, as shown in Fig. [IJa). This assumption allows us to decouple Eq. ^ into 
two separate second-order differential equations for the amplitudes uq = f{z)exp{—ikfz) and vq = g{z) exp(—ikfz), 
by neglecting terms proportional to 9^/ compared to terms in kfdzf, and the same for g. Each of the resulting 
equations has two independent solutions of the form exp[±ifcz], where k ^ [v'jieo — hkfvY — Aq{v^ — w^)]^/^ with 
Vf = Hkf /m being the Fermi velocity, cq being the energy of the Andreev bound state [the dashed line in the Fig. [2] 
inset], and Aq being the bulk order parameter. The expression under the square root is zero at the pair-breaking 
velocity v w Vpb ~ Vf\Ao\/2fi in the BCS regime. For smaller velocities k is imaginary and hence the Andreev bound 
state is localized. The size of the bound state is thus It ~ ^/\k\, which diverges at the pair-breaking velocity. Our 
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FIG. 2: (a) Profile of the modulus of the order parameter |A(x)j for l/kf-a — —0.5 and v — (solid blue curve), 0.2c (solid green 
curve with cicles), 0.3c (dashed red curve) and 0.38c (dot-dashed black curve), obtained from the time-independent calculations 
of a homogeneous superfluid. (b) Corresponding plots of the lowest Andreev bound state (specifically [uo(a;)|^). Inset: the 
spectrum of energy levels in a homogeneous gas containing a soliton. The solid line shows the excitation gap in the continuous 
spectrum of the extended Bogoliubov states in a homogeneous Fermi gas as a function of velocity for l/fc/ffl = —0.5. The 
shaded area above represents the continuum spectrum for the homogeneous infinite system. The dashed line shows the energy 
of the lowest energy level, which is the Andreev bound state shown in (b). 

calculations suggest that the soliton solution ceases to exist at a velocity slightly below Vpb where the excitation gap 
is still finite. 

In Fig. [3ja) we plot the energy of the soliton Eg as a function of v (the energy dispersion) for various values of 
l/kja. In the BEC limit, the energy is given by the well-known Gross-Pitaevskii result [Tsl. [l9l [35l| 

i?,oc(l-«Vc')", (3) 

where a — 3/2. For 1/kfa = 1 (black dotted curve), our numerical data closely follows this result. At unitarity (red 
dashed curve), the energy initially decreases smoothly with increasing v, as for 1/kfa = 1, but then reaches a local 
minimum, subsequently briefly increases, before decreasing again rapidly. Finally the energy dispersion truncates at 
a small but finite value of Eg [36j . In a previous publication, some of us predicted that at unitarity a = 2 in Eq. 
(|31) [Tfij . This is approximately true for small v. The presence of a local minimum in the energy dispersion is a crucial 
point, as we now explain. 

In general, the soliton energy Eg is a function of /i and [l5|. In Fig[IJc), the soliton begins at a distance Xq 
from the trap center with v = 0. The soliton may move closer to the center of the trap, to larger /.t, maintaining a 
constant Es, by increasing u^. However, once the soliton reaches the local minimum in the energy dispersion, it can 
no longer move to larger /i by changing w^. Hence, in a real experiment, the soliton must decay when v reaches the 
position of the local minimum, not when it reaches the point where the curve truncates. 

Following this line of reasoning, we may reach another important conclusion. If the energy of the soliton were to 
decrease monotonically to zero with increasing \v\, as in the Gross-Pitaevskii equation, decay of the kind shown in 
Fig- HJc) would be impossible. This is because the soliton could always maintain a constant E'g, however much /i 
increased, by increasing v closer to the maximum v where Eg goes to zero. Hence, the existence of the decay process 
shown in Fig. [TJc) implies a truncation in the energy dispersion, or at least a non-monotonic energy dispersion. 

When we reduce 1/kfa from to —0.2 [green solid curve in Fig.jSja)], we again find a local minimum in the energy 
dispersion, but the curve truncates at a much larger value of Eg. We now obtain a ~ 2.2 in Eq. The minimum 
value of Eg is even larger for 1/kfa — —0.5 (blue dot-dashed curve), and a has increased to 2.8, but now the energy 
dispersion does not contain a minimum. Instead, we see a broad bump in the curve as the soliton approaches the 
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FIG. 3: a) Energy of soliton Es in units of the Fermi energy Ef as a function of soliton velocity v for 1/fc/a = —0.5 (blue dot- 
dashed curve), —0.2 (green solid curve), (red dashed curve) and 1 (black dotted curve), obtained from the time-independent 
calculations of a homogeneous superfluid. The curves of numerical results for 1/fc/a = —0.5, —0.2 and are fitted to thin solid 
curves of the form Es oi [l — u^/c^)". The solid horizontal lines indicate Vpb for 1/fc/o = —0.5, —0.2 and from left to right, 
(b) Corresponding plot for the phase jump across the soliton A(j>- 

maximum velocity. 

The non-monotonic slope of the energy dispersion in the BCS regime can be tentatively explained by considering 
the Friedel oscillations. These arise due to scattering of fermions at the Fermi energy by the soliton, and appear 
as ripples in the density or order parameter profile either side of the central minimum of the soliton [as shown in 
Fig. [TJa), for example]. For Friedel oscillations to appear, the atomic wave length at the Fermi energy, A w 27r/fc/, 
should be larger than the size of the soliton. When the velocity of the soliton increases, its size grows and can exceed 
A, causing the loss of the Friedel oscillations. For example, we find that the size of the bound state 4 and A are 
comparable at v ^ 0.3c for 1/kfa — —0.5. For v — 0, 0.2c and 0.3c, the Friedel oscillations are almost identical in 
the profile of the order parameter [blue solid curve, green solid curve with circles and red dashed curve in Fig.[21[a), 
respectively], although the depth of the soliton varies dramatically. In contrast, if v increases to 0.38c, the Friedel 
oscillations disappear and the order parameter acquires a smooth broad profile [black dot-dashed curve in Fig.[2Ua)]. 
Due to the disappearance of the Friedel oscillations, we would expect different behaviour in the energy dispersion, as 
observed in Fig. [3]Ja). Notice that the existance of a small range of v where Eg increases with v implies that, if a 
soliton were created by phase imprinting in this range of v, the snake instability would be absent f37| . 

Figure [3fb) shows the equivalent plot to Fig. [Jl^a) for the phase jump across the soliton A<j). In the BEC limit, the 
Gross-Pitaevskii equation predicts that A(/) — —2 arccos {v/c). For 1/kfa — 1 (black dotted curve), our numerical data 
again agrees well with the Gross-Pitaevskii prediction. At unitarity (red dashed curve), A0 varies linearly with v for 
small V, as some of us predicted analytically in a previous publication [l6j . However, as v approaches the maximum 
velocity we see a departure from this linear behaviour: in a similar manner to the energy dispersion, the A(f)(v) curve 
forms a local minimum, then briefly rises before decreasing rapidly and terminating. This plot makes the important 
point that the soliton still has a finite phase jump when the energy dispersion terminates. This is expected because, 
as shown in Fig. [2Ua), the energy dispersion truncates at finite energy. In the case of a soliton moving with a nonzero 
velocity v, a finite phase step implies finite energy by virtue of the kinetic energy associated with a superfluid current. 
When we decrease 1/fc/a from to —0.2 we find a similar variation of A0 with v (green solid curve), except that the 
curve now truncates at a larger A0. As for the energy dispersion, at 1/fc/a = —0.5 (blue dot-dashed curve) we find a 
broad bump instead of the local minimum. 

Using the above results we plot the maximum soliton velocity Vm in the homogeneous gas as a function of 1/kfa 
as red squares in Fig. [J] As explained earlier, in a real or numerical experiment, we predict that the measured Vm is 
given by the position of the local minimum in the energy dispersion, not by the the point where the curve truncates, 
and hence we define Vm to be the position of that minimum. For 1/kfa — —0.5 there is no local minimum in the 
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FIG. 4: Dashed curve: sound speed c in a homogeneous superfluid as a function of l/kfU. Solid curve: bulk pair-breaking 
velocity Vpb- Red squares: maximum soliton speed Vpb obtained from the time- independent calculations of a homogeneous 
superfluid. Green crosses (plus sign): maximum soliton speed obtained from time-dependent calculations of a superfluid in a 
trap with frequency ujx = 2-k x 50 {2-k x 25) rad s^^. Blue downward (upward) pointing triangles: soliton speeds before (after) 
the collisions presented in the middle row of Fig. [5l 

energy dispersion, but we take Vm to be the velocity at the "broad bump" for consistency. We compare this data 
to the bulk sound speed c (dashed curve) and bulk pair-breaking velocity Vpb (solid curve) for a homogeneous gas, 
given by mv'^^^ — •\/A^+/? — /x. For 1/fcya — 1, the maximum soliton velocity is given by c, as expected from the 
Gross-Pitaevskii equation fl^. The pair-breaking speed is much larger and hence not relevant because the pairs are 
strongly bound as molecules. However, for l/kfu < 0, the pair-breaking speed determines the Landau critical velocity, 
and hence limits Vm- Thus, at unitarity we find a Vm which is slightly smaller than Vpb- As we decrease ^/kfa we find 
that Vm closely follows the Vpb curve. In fact, the data points lie slightly below Vpb, and the deviation becomes larger 
as l/kfa decreases. 

On the same graph we also include a number of points (black crosses) showing the critical velocity predicted by 
the time-dependent simulations of a trapped superfluid. The points lie slightly below the predictions of the time- 
independent calculations, the relative deviation being less at unitarity than at smaller values of 1/fc/a. This small 
deviation is a finite-size effect, which probably occurs because the pair size for small l/kfu becomes comparable to the 
size of the cloud. Hence it would be exceedingly difficult to probe soliton dynamics near the bulk pair-breaking speed 
in a real experiment. We include one point (green plus sign) in Fig. 2] for oj^ = 2Tr x 25 rad s~^ and l/kfa = —0.5 to 
illustrate that the critical velocity in the time-dependent simulations approaches the prediction of the time- independent 
calculations as the trap frequency is reduced. 

IV. SOLITON COLLISIONS 

We place two black solitons in a trapped Fermi superfluid, with a displacement of ±Xo from the trap centre, and 
then evolve in time. The solitons are accelerated by the harmonic potential and collide at the trap centre. We may 
increase the speed of the collision by increasing Xq. Figure [5] shows the evolution of the density proflle with time for 
different Xq and l/kja, with the other parameters as in Fig. [TJ In the upper (middle) row Xq = 5.1 (2.1) /im, and 
l/kfa decreases from left to right (see caption). 

Let us begin with the BEG regime l/kfa = 1, shown in the left-hand column. For Xq = 5.1 fj,m [Fig. [5ja)], the 
two solitons pass through each other without changing their form, and subsequently slow down and come to rest at 
X ~ ±Xo, indicating that no energy has been lost by the solitons and hence the collision is elastic. For Xq = 2.1 /im 
[Fig.[5l^f)], the solitons behave differently. In this case, the solitons slow down as they approach one another, come to 
a halt, and then move away without crossing. This is known as a black collision, whereas a crossing of two solitons, as 
shown in Fig.[5Ua), is known as a grey collision. We again note that the two solitons return to their original positions 
after the collision, indicating that the collision is elastic. This behaviour of grey collisions for high impact speeds, and 
black collisions for low impact speeds, has already been observed in the CP equation [SQ] . 
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FIG. 5: Collisions of solitons across the BEC-BCS crossover. The top (middle) row shows grey-scale plots of the density profile 
n{x,t) (black=high) in a superfluid, in which two solitons begin at rest with a displacement Xo = ±5.1 (2.1) fj,m from 
the trap centre, cu^ = 2tv x 50 rad s~^, L± = 3.3 ^m and the peak density Up = 1.8 x 10^* m~'^. The bottom row shows 
the evolution of the lowest Andreev state (specifically |uo(a;,t)|^) in grey-scale (black=high) during the collisions shown in the 
middle row. 1/fc/ffl = 1.0, 0.2, 0.0, —0.35 and —0.5 in each column from left to right. 

When we reduce l/kfu to 0.2, the behaviour is similar. We again observe a grey collision for Xq = 5.1 /im [Fig.[5t^b)], 
and a black collision for Xq = 2.1 /xm [Fig. [5l^g)]. However, on close inspection we observe that, after the collision, the 
solitons come to rest at a value of > Xq, indicating that the solitons have a greater speed after the collision than 
before. This is particularly noticeable for the slower collision [Fig. [Sl^g)]. Since, counter-intuitively, fast solitons have 
less energy than slow solitons (as shown in Fig [3]), we conclude that the solitons have lost energy due to an inelastic 
collision. We also observe small ripples in the density which emanate from the collision, showing that the energy lost 
by the solitons is converted into sound. 

We interpret this effect as being due to the evolution of fermionic quasiparticles localised in the solitons. In the 
bottom row of Fig. [5l we plot the evolution of the lowest Andreev state (specifically \vo{x,t)\'^) during the slow 
soliton collisions shown in the middle row of Fig. [5j As the solitons approach, the Andreev states see a double-well 
potential, and subsequently a deep single- well (if the collision is grey), causing them to oscillate and breathe in a 
non-adiabatic way. In the BEC regime this effect is negligible, because the Andreev states make a small contribution 
to the overall density (/ uq dx dy dz ^ J vq dx dy dz, for example), and because the solitons repel each other at slow 
speeds. Consequently the slow collision causes no disruption to the lowest Andreev state for 1/fc/a — 1.0 [Fig. [^k)]. 
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However, for 1/kfa = 0.2 we observe oscillations in the amplitude of the lowest Andreev state [Fig. EJl)]. At one 
instant, shortly after the collision, |wo(a;,t)|^ reduces almost to zero. This means that particles are being transferred 
between different eigenstates of the Bogoliubov spectrum, with a coupling between the localized Andreev state and 
the states in the continuum. These transitions are associated with density and phase oscillations which eventually 
cause an emission of sound and loss of energy from the soliton. This effect is enhanced for low impact-collisions, 
because the fermionic quasiparticles have more time to move in response to the change in potential. The key role 
played by fermionic quasiparticles in this process is also confirmed by the fact that the same collisions are found to be 
elastic when a purely bosonic density functional theory (i.e., a nonlinear Schrodinger equation for the pairing field) is 



used in place of the TDBdG equations [32 1 



As we further reduce 1/kfa, the collisions become increasingly inelastic. At unitarity [Figs. [S{c) and (h)], the 
collisions create a great deal of sound and the solitons subsequently come to rest at much larger than ATq, particularly 
for Xq ~ 2.1 /im [Fig. (S^h)]. The inelastic collision is associated with very pronounced oscillations in [uol^^iOl 
[Fig. ETrn)]. We also observe that the collision for ATo = 2.1 ^m is no longer black, due to the increasing mass of the 
solitons [l3- When 1/kfa reaches —0.35, the sohtons are destroyed by the collision for — 5.1 ^m [Fig.[5jd)], but 
survive if Aq = 2.1 fim [Fig.lSfi)]. Solitons can be destroyed by the inelastic collisions if their residual energy after the 
collision is less than the minimum soliton energy [38j . For l/kfU = —0.5 and Xq = 5.1 /zm, the minimum energy (or 
maximum speed) is reached before the solitons even collide [Fig.[SJe)]. For Xq = 2.1 fj,m, the solitons collide but then 
are immediately destroyed [Fig. [5l^j)]. Consequently the lowest Andreev state is no longer localised after the collision 
[Fig. Ho)]. 

We quantify how inelastically the solitons collide in Fig. 01 Here we plot the soliton speed immediately before and 
after the soliton collisions presented in the middle row of Fig. [5] as pairs of downward and upward-pointing triangles 
respectively. For 1/kfa = 1, the collision is elastic, and hence the downward and upward-pointing triangles lie on top 
of each other. For 1/kfa = 0.2, the triangles are a small distance apart, indicating a slightly inelastic collision. This 
distance increases as 1/kfa is reduced to 0, showing that the collision is becoming increasingly inelastic. However, 
the speed of the soliton after the collision is still far below the maximum speed observed for a single soliton, shown by 
the black crosses. However, when 1/kfa is reduced to —0.35 the final soliton speed is very close to the critical value. 
This explains why collisions for larger Xq or smaller 1/kfa destroy the solitons. 



V. CONCLUSIONS 



We have shown that in a two-component superfluid Fermi gas soliton solutions cease to exist as the soliton ap- 
proaches the bulk pair-breaking velocity Vpb on the BCS side of the resonance. At this point, the variation of soliton 
energy with velocity (the energy dispersion) truncates at a finite value of energy and phase jump across the soliton. 
Close to Vpb, the energy dispersion shows interesting structure, and sometimes a local minimum, which can reduce 
the maximum soliton velocity as measured in a real or numerical experiment. The presence of a harmonic trap causes 
a further reduction in the maximum observed soliton velocity. If a soliton is accelerated to its maximum velocity by 
an external potential, it abruptly disappears and its energy is converted into sound. Solitons may also be accelerated 
by inelastic collisions with other solitons, and the energy lost by the solitons is emitted as sound. This can destroy 
the solitons if the energy lost as sound reduces the soliton energy below the truncation point in the energy dispersion. 
We find that soliton collisions become increasingly inelastic as we move from the BEC to BCS regime. 

On the one hand, our results impose some limitations on experiments aimed to observe soliton oscillations and 
collisions in Fermi gases. Solitons on the BCS side of the resonance must be prepared close to the trap centre and 
accelerated gently in order to avoid their decay into sound. On the other hand, it is encouraging that solitons in 
fermionic superfluids behave differently to those in bosonic superfluids, and that this physics is not captured by 
bosonic models. The solitons in fermionic superfluids are sensitive to the fermionic degrees of freedom and so may be 
used as a tool to further characterize these gases. 
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